Theory of correlations in strongly interacting fluids of two-dimensional dipolar bosons 
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Ground-state properties of a two-dimensional fluid of bosons with repulsive dipole-dipole inter- 
actions are studied by means of the Euler-Lagrange hypernetted-chain approximation. We present 
a self-consistent semi-analytical theory of the pair distribution function g(r) and ground-state en- 
ergy of this system. Our approach is based on the solution of a zero-energy scattering Schrodinger 
equation for the "pair amplitude" \J g(r) with an effective potential from Jastrow-Feenberg cor- 
relations. We find excellent agreement with quantum Monte Carlo results over a wide range of 
coupling strength, nearly up to the critical coupling for the liquid-to-crystal quantum phase transi- 
tion. We also calculate the one-body density matrix and related quantities, such as the momentum 
distribution function and the condensate fraction. 
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I. INTRODUCTION 

Although experiments on ultra-cold atomic gases are 
performed in an extremely dilute regime, inter-particle 
interactions can play an important role in determining 
their behavior £Q. In the ultra-cold regime the van der 
Waals interaction between atoms can be replaced by a 
Fermi pseudo-potential whose sign and strength can be 
tuned at will by means of a Feshbach resonance pQ . 

Recent experimental advances (2 in cooling down 
atoms with permanent dipole moments and polar 
molecules have made it possible to study quantum gases 
with long-range interactions. These systems have at- 
tracted a large fraction of theoretical and experimental 
interest 0] . The fact that the strength and sign of the 
short range part of the inter-molecular interaction can be 
fully controlled by means of Feshbach resonances makes 
polar molecules an excellent system for studying the in- 
terplay between short-range and long-range interactions. 
Moreover the short-range component of the scattering 
potential can be completely switched off, leaving us with 
a purely dipolar system. The dipole-dipole interaction 
is peculiar in the sense that it is both long range and 
anisotropic. The latter property is interesting in par- 
ticular since the dipolar interaction can become either 
attractive or repulsive depending on the geometry of the 
external confining potential. 

Magnetic dipole moments of atoms lead to an 
anisotropic and long-ranged dipole-dipole interaction and 
experiments have shown that the dipole-dipole interac- 
tion can influence the shape and stability of a quantum 
gas [5]. Furthermore, there has been a lot of experimental 
progress in stabilizing and cooling dipolar molecules [6] 
with an electric dipole moment which is larger than their 
magnetic dipole moment, making the dipolar interac- 
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tion the dominant one. Dipole-dipole interactions be- 
tween molecules can be thus controlled by external elec- 
tric fields. 

A number of theoretical and computational stud- 
ies have addressed ground-state properties of two- 
dimensional (2D) repulsive dipolar bosons at zero and 
low temperatures [7HS]- Quantum Monte Carlo (QMC) 
studies [5] at zero temperature have covered, in particu- 
lar, the whole range of coupling strength up to the crys- 
tallization transition. 

In this Article we present a theoretical study of ground- 
state properties of a polarized dipolar fluid of bosons 
with average density n, tightly trapped in the direc- 
tion of polarization. The system is effectively 2D and 
the dipole-dipole interaction is isotropic. The main fo- 
cus of this Article is on the so-called pair distribution 
function (PDF) [10] g(r), which is defined so that the 
quantity 2i:ng(r)dr gives the average number of dipoles 
lying within a circular shell of radius dr centered on a 
"reference" dipole sitting at the origin. We present a 
self-consistent semi-analytic theory of the PDF, which 
incorporates many-body correlation effects, allowing us 
to explore the physics of the system at strong coupling. 
Our approach, which is based on the so-called Euler- 
Lagrange hypernetted-chain (HNC) approximation [TT1 - 
IT5] . involves the solution of a zero-energy scattering 
Schrodinger equation with an effective potential [TOll8| . 
Following Ref. [T7J, we transcend the HNC approxima- 
tion by taking into account higher-order {i.e. three-body) 
correlations in an effective way. We also calculate, within 
the lowest-order HNC formalism, the one-body density 
matrix, the momentum distribution function, and the 
condensation fraction. 

A similar study has recently been carried out by Huf- 
nagl et al. [12]. The authors of this work have taken 
into account both the anisotropy of the dipolar interac- 
tion and the short-range part of the repulsive potential, 
but have neglected three-body correlations, thereby ob- 
taining very good agreement with QMC data only in the 
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weak-coupling regime. 

The outline of this paper is as follows. In Sec. [TT] we 
briefly review the main properties of the dipole-dipole 



interaction. In Sec. Ill we present a formally exact zero- 
energy scattering equation for the pair amplitude \J g(r) 
and introduce the approximations that we employ for the 
evaluation of the effective potential. Our main numerical 
results are presented in Sec. |IV| Finally, Sec. |V] summa- 
rizes our main conclusions. 



II. DIPOLAR INTERACTIONS 

The interaction energy between two particles with 
identical dipole moments aligned along the unit vectors 
rii and n 2 and placed at a distance r\ 2 from each other 
is given by the following well-known expression: 



Vdd(»"i 



Cdd 

47T 



(ni • n 2 ) 



3(ni • r 12 )(n 2 ■ r 12 ) 



(1) 

Here Cdd is the dipole-dipole coupling constant, which 
depends on the microscopic origin of the interaction: e.g., 
it is d 2 /eo for particles with permanent electric dipole d 
and [IqM 2 for particles with permanent magnetic dipole 
M (here eo and fJ-o a re the permittivity and permeability 
of vacuum, respectively). An electron- hole bilayer in the 
exciton condensate phase [5D] is another realization of 
a system with dipolar interactions. In this case Cdd = 
e 2 d 2 /e, where — e is the electron's charge, d is the inter- 
layer separation, and e is the dielectric constant of the 
host semiconductor. 

For the polarized case, where all dipoles align in the 
same direction, Eq. ([I]) simplifies to 



^dd(?"l2) = 



Cdd l-3cos 2 (6») 
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This interaction has two important features. It is long- 
ranged, i.e. it decays like 1/r 3 at large distances, and 
furthermore, it is anisotropic. In particular, it is repul- 
sive for dipoles aligned side-by-side (8 = tt/2) and is at- 
tractive for dipoles aligned head-to-toe (9 = 0). In this 
work we study 2D systems of polarized dipoles, so that 
the inter-particle interaction is isotropic: 
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(3) 



III. THEORY 



We consider a 2D fluid of N point dipoles with bosonic 
statistics. The first-quantized Hamiltonian can be writ- 
ten as [5]: 



H 



2m 



Cdd 
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(4) 



m being the mass of a nano-particle dipole. The ground- 
state properties of the Hamiltonian Q are governed by 
a single dimensionless parameter: 



nr. 



! 
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where n is the average density and r = mC^dl (4irh 2 ) 
is a characteristic length scale, which is typically of the 
order of a few Angstroms. 

In order to calculate the ground-state properties of the 
Hamiltonian Q , we use the HNC [TTHT5] approximation 
at zero temperature. In what follows we first present our 
theory at the simplest level (which works well in the limit 
7 <C 1) and then transcend it to obtain accurate results 
at strong coupling (7 3> 1). 

With the zero of energy taken at the chemical poten- 
tial, the formally exact differential equation for the pair- 
correlation function can be written as 1151 1161 
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The effective scattering potential is 

V eS {r)=v dd (r)+W B (r) 



(6) 
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where i>dd(f) is the bare repulsive dipole-dipole interac- 
tion in Eq. ^ and W B {r) is defined, at the level of the 
so called "HNC/0" approximation, by the following equa- 
tion |2T] : 



W B (k) 



e{k) 



[2S(k) + 1] 



S(k) - 1 
S(k) 



(8) 



In writing Eq. d8| we have introduced the Fourier trans- 
form (FT) W B (W) of W B (r) according to 



FT[F(r)] = n J d 2 r F{ r) exp (ik ■ r) 



(9) 



Furthermore, e(k) — h 2 k 2 /(2m) is the single-particle en- 
ergy and S(k) is the instantaneous or "static" structure 
factor 10 . S(k) = 1 + FT[g(r) - 1]. Eqs. form a 

closed set of equations, which can be solved numerically 
self-consistently to the desired degree of accuracy. 

When 7 ~ 1 the simplest approximation for W B (r) 
in Eq. ^ is inadequate. Improvements on Eq. ^ for 
a Bose fluid can be sought in two directions |22| . The 
HNC/0 may be transcended by the inclusion of i) low- 
order "elementary" diagrams and ii) three-body Jastrow- 
Feenberg correlations. 

The contribution from three-body correlations is [22] 



Wi 3 \k) 



JL_ r d 2 q 
An J (2tt) 2 

[E(p) + E(q)]u 3 (q,p,k)} 



S(p)S(q)u 3 (q,p,k){v 3 (q,p,k) 

(10) 



In Eq. ([10]) p = -(q + k), E(k) = e(k)/S{k) is the 
Feynman-Bijl excitation spectrum (see Sect. IV), 



v 3 (q,p,k) = (h 2 /m)[k-px(p)+k-q X (q)+p-qx(q)} , (H) 
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and 



u 3 (q,P,k) = - 



(h 2 /2m) 
E(k) + E(p) H 

[k ■ px(k)x(p) + 
k ■ qx(k)x(q)} ■ 



E(q) 

p ■ qx(p)x(q) 



(12) 



In Eqs. ([n}-(|l2]) X (k) = 1 - 1/S(k). 

We have taken into account higher-order terms that 
are missed by the HNC/0 approximation by assuming 
that they lead to corrections to the scattering potential 
V e g(r). Using the theory developed by Apaja et al. JS3] 
we have supplemented Wb(&) in Eq. (Tsl) by the inclusion 
of the three-body potential W B (k): 

W B (k) -> W B {k) + a( 1 )W B 3) (k) . (13) 

If a(-y) is set to unity, the r.h.s. of Eq. ( [T3| defines the so- 
called "HNC/3" approximation. It has been shown [T7] 
that higher-order corrections beyond HNC/3 can be ef- 
fectively taken into account by introducing a weighting 
function a(-f) > 1. This approximation will be termed 
"HNC/a3". 

The functional dependence of the weighting factor a(j) 
on the coupling constant 7 can be fixed, for example |17j , 
by requiring that the ground-state energy per particle ex- 
tracted from the PDF calculated within the HNC/a3 ap- 
proximation matches exactly the corresponding quantity 
calculated by QMC simulations [8]. Following this pro- 
cedure we have calculated numerically the function 0(7) 
and produced a convenient analytical parametrization of 
it in the interval 1 < 7 < 256: 



"(7) 



3.26 exp (-O.267 



0.56\ 



(14) 



Numerical results for the PDF calculated within the 
three approximations described in this Section (HNC/0, 
HNC/3, and HNC/a3) will be illustrated in Sect. [TV 
and severely tested against accurate QMC data by As^ 
trakharchik et al. [5]. Practical recipes on how to solve 
Eq. (JsJ) are discussed in detail in Ref. [13] . 



A. The one-body density matrix 

Although the main focus of this Article is on the PDF, 
in this Section we briefly discuss how to calculate the 
one-body density matrix (1BDM) and other interesting 
quantities on the basis of the simplest scheme discussed 
earlier, i.e. the HNC/0. 

Once the PDF g(r) is obtained from the formalism out- 
lined in the previous Section, one can use it to find the 
1BDM and hence the condensate fraction and the mo- 
mentum distribution function. The 1BDM of a homoge- 
nous quantum liquid is related to its many-body wave 
function \&(t"i, r%, ■ ■ ■ , t*at) by the following relation: 

p(r) = J d 2 r 2 ■ ■ .d 2 r N ^*(r,r 2 , ■ . ■ ,r N ) 

x 9(0,r 2 ,...,r N ) , (15) 




FIG. 1: (Color online) The effective scattering potential 
V e ff(r) of a 2D fluid of dipolar bosons [in units of 'yti 2 / (mr^)] 
as a function of n 1 / 2 ?-. Different curves correspond to differ- 
ent values of the dimensionless coupling constant 7. All the 
curves in this plot have been calculated within the HNC/a3 
approximation. 



where we have assumed the following normalization 
/ d 2 ri . . .d 2 r N \^(r 1 ,r 2 , . . . ,t"at)| 2 = N, N being the to- 
tal number of particles. 

In a homogeneous system the value at the origin of 
the 1BDM gives the particle density, p(0) = n, while its 
long-distance behavior gives the condensate fraction no 
(i.e. the fraction of particles occupying macroscopically 
the zero-momentum state, here measured in units of the 
total density n), uq = p[r — > 00) /n. Furthermore, the 
momentum distribution function n(k) is also related to 
the 1BDM via Fourier transformation: 

n(k) = nn (2Tr) 2 6{k) + FT[p(r)/n - n Q ] . (16) 

Within the HNC/0 formalism, the 1BDM reads [HIM] 



p(i 



nnoe 



JV m (r) 



(17) 



where the N ww (r) nodal function is given (in Fourier 
transform) by 

^ww(fc) = [S wd (k) - 1] [S wd (k) - 1 - AU(£:)] ■ (18) 

Here 5 w d(fc) and iV wd (fc) could be obtained from the so- 
lution of the following coupled equations 



and 



Here 



N wd (k) = [S wd (k) - 1] [S(k) - 1 - N(k)} 
S wd {k) = l + FT[g wd (r)-l] , 

<?wd(r) = f(r)e N «^ . 
[S(k) l] 2 



N(k) 



S(k) 



(19) 
(20) 

(21) 
(22) 
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FIG. 2: The pair distribution function g(r) of a 2D fluid of dipolar bosons as a function of n 1 ^ 2 r for different values of 7. Dashed, 
dotted, and solid lines represent theoretical results based on the HNC/O, HNC/3, and HNC/a3 approximations, respectively. 
Filled circles denote QMC results from Ref. [8]. 



is the conventional nodal function and the correlation 
function f(r) is given by f(r) — g(r) exp [— N(r)]. 

Once self-consistency in the solution of Eqs. (fl9|-(21) 



is achieved, we can find the condensation fraction via 

n = exp (2i? w - R d ) , (23) 



where 



n / d 2 r[g wd (r) - 1 - N wd (r)] 



d 2 r[g wd (r) - l]iV wd (r) , 



(24) 
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and 



<fr \g{r) - 1 - N(r)] 
d 2 r [g(r) - 1] N(r) . 



(25) 



All the necessary ingredients to find the 1BDM and the 
momentum distribution function are now available. 



IV. NUMERICAL RESULTS 

We now turn to a presentation of our main numerical 
results. 

Fig. [I] illustrates the dependence of the effective scat- 
tering potential V c g(r) calculated at the HNC/a3 level 
on the interaction strength 7. Note that at strong cou- 
pling the behavior of V c ff(r) is highly oscillatory indicat- 
ing the emergence of short-range order with increasing 7 
and the incipient quantum phase transition to an ordered 
crystalline phase at strong enough coupling. According 
to the QMC study by Astrakharchik et al. [5], 7 ~ 290 is 
the critical value at which this quantum phase transition 
occurs. 

The degree of agreement between our theoretical re- 
sults for the PDF and corresponding QMC data is illus- 
trated in Fig. [2] The inclusion of effectively-enhanced 
three-body correlations (HNC/a3 theory) yields an ex- 
cellent agreement with QMC data over the entire range of 
coupling constants where the fluid phase is stable. Note 
that when two point dipoles are in close proximity to 
each other, the PDF is determined by the solution of a 
two-particle Schrodinger equation in the r — > limit. Ac- 
cording to Kimball's cusp condition [25] one indeed finds 
g(r -> 0) ~ v /r/r exp(-4 v /r /r). 

The ground-state energy per particle £qs of the sys- 
tem can be calculated by means of the well-known 
integration-over-the-coupling-constant algorithm [10J: 



n 

2 jo 
nC dd 



dX / d 2 r v dd (r)g x (r) 



dX 



00 dri^ 



(26) 



where g\{r) is the PDF of an auxiliary system 
with dipole-dipole interactions of the form £'(r) = 
Xv dd (r) — XC dd /(4irr 3 ). In practice the integration over 
A is carried out by integrating over ^7. In Fig. [3] we 
have reported the ground-state energy calculated at the 
HNC/0 level in comparison with the QMC result [5]. (By 
construction, the ground-state energy calculated at the 
HNC/a3 level coincides with the QMC result.) 

Figs. Hand [5] illustrate the static structure factor S(k) 
and the (approximate) excitation spectrum E{k) for sev- 
eral values of 7, respectively. As the coupling constant 
is increased, correlations get stronger and the height of 
the first-neighbor peak in S(k) increases. Note that S(k) 
vanishes linearly in the long- wavelength limit. 




FIG. 3: (Color online) The ground-state energy (per particle) 
of a 2D fluid of dipolar bosons £gs in units of yti 2 / {mr^) as a 
function of interaction strength 7. The dashed line represents 
our HNC/0 results. The solid curve is the parametrization 
formula of the QMC data reported by Astrakharchik et al. [8]. 
We have not plotted our HNC/a3 results since they coincide 
by construction with the QMC results (we remind the reader 
that enforcing this constraint allows us to find the weight- 
ing factor 0(7) that enhances triplet correlations at strong 
coupling) . 



We calculate the excitation spectrum E(k) from the 
approximate Feynman-Bijl relation |26j 



E(k) 



e{k) 
S(k) 



(27) 



which is believed to be asymptotically exact for k — > 
and provides a rigorous upper bound on the excitation 
energy at finite k (since it represents the average en- 
ergy of the excitations which couple to the ground state 
through the density and thus necessarily exceeds the min- 
imum excitation energy [27]). Since S(k — > 0) oc k, one 
finds collective excitations with an acoustic dispersion at 
small k. In agreement with the QMC data [5], we observe 
that deviations from the linear phononic behavior at 
small k start very soon and that a strong roton minimum 
appears at finite k as the interaction strength increases. 
The corresponding roton gap decreases upon increasing 
the interaction strength. The structure of the phonon- 
roton peak of 2D dipolar bosons has been recently stud- 
ied by Mazzanti et al. |28j by employing the dynamical 
structure factor S(k,u>). Moreover, multi-particle excita- 
tions, which are absent in the Feynman-Bijl theory, have 
also been taken into account by the same authors 28J. 

Finally, in Figs.[6]|8]we illustrate HNC/0 results for the 
1BDM p(r), the momentum distribution function n(fc), 
and the condensate fraction uq, respectively. 

Our HNC/0 results for the 1BDM are in good agree- 
ment with the corresponding QMC data [29 at large r 
(see Fig. [6| , but deviate considerably from them at short 
distances. Specifically, we do not recover the exact result 
p(0) = n. This is a well known deficiency of the HNC/0 
approximation [24 . 
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FIG. 4: (Color online) The instantaneous structure factor 
S(k) of a 2D fluid of dipolar bosons as a function of k/n 1 / 2 
for various values of 7. Data in this plot refer to the HNC/a3 
approximation. 



FIG. 6: (Color online) The one-body density matrix p(r) of 
a 2D fluid of dipolar bosons (in units of n) as a function of 
n 1 ^ 2 r for various values of 7. All the theoretical results in 
this plot (lines) refer to the HNC/0 approximation. Symbols 
refer to QMC results from Ref. [SS]. 





FIG. 5: (Color online) The upper bound E(k) on the ex- 
citation spectrum of a 2D fluid of dipolar bosons in units of 
7/i 2 / (mr 2 ) is plotted as a function of fc/n 1//2 for various values 
of 7. The (black) dashed line is the noninteracting parabolic 
spectrum e{k) = h 2 k 2 /(2m). All the other lines refer to the 
HNC/q3 approximation. Symbols indicate the prediction for 
E(k) obtained from the QMC results [8] for S{k). 



The condensate fraction of 2D fluids of bosons is an 
intriguing quantity. It is well known [30 that 2D bosons 
with ln(r) interactions, despite displaying superfluidity, 
are characterized by a vanishing condensate fraction (at 
all temperatures). On the other hand, 2D bosons with 
1 jr interactions display a finite value of n at zero tem- 
perature [3T]. Fig.[8jillustrates the condensate fraction of 
2D bosons with l/P interactions, which is finite and in 
very good agreement with that predicted by QMC sim- 
ulations [H], even at strong coupling. In the limit 7 — !• 
one expects [35] n m 1 — 1/| I117I. 



FIG. 7: (Color online) The momentum distribution function 
n(k) of a 2D fluid of dipolar bosons as a function of k/n 1 ^ 2 
for various values of 7. 

V. SUMMARY AND CONCLUSIONS 

In summary, we have presented a self-consistent semi- 
analytical theory of correlations in strongly interacting 
fluids of two-dimensional dipolar bosons. Treating in 
an effective manner high-order correlations beyond the 
so-called HNC/0 approximation [21], we find excellent 
agreement between our results for the pair distribution 
function and quantum Monte Carlo data [8]. Sticking to 
the simplest HNC/0 method, we also calculate the one- 
body density matrix, the momentum distribution func- 
tion, and the condensate fraction. 

Our results are also extremely useful to study two- 
dimensional dipolar fermions, for which quantum Monte 
Carlo data have recently appeared in the literature |33j . 
Fermions can be easily treated within our formalism by 
including a statistical contribution to the effective scat- 
tering potential jTS] (which is crucial to enforce the Pauli 
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FIG. 8: (Color online) The condensate fraction no of a 2D 
fluid of dipolar bosons as a function of 7. Empty circles repre- 
sent the HNC / results (the line connecting the empty circles 
is just a guide to the eye). Filled circles are QMC data from 
Ref. 0. 



principle on the on-top value of the pair distribution func- 
tion for fermions with antiparallel spin). 

A straightforward generalization of our theory is also 
useful to study density instabilities in two-dimensional 
dipolar fermions going beyond the Singwi-Tosi-Land- 
Sjolander approximation |34j . 
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